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ABSTRACT 


In  this  paper  we  present  a  high  order  weighted  essentially  non-oscillatory  (WENO)  scheme 
for  solving  a  multi-class  extension  of  the  Lighthill-Whitham-Richards  (LWR)  model.  We 
first  review  the  multi-class  LWR  model  and  present  some  of  its  analytical  properties.  We 
then  present  the  WENO  schemes,  which  were  originally  designed  for  computational  fluid 
dynamics  problems  and  for  solving  hyperbolic  conservation  laws  in  general,  and  demonstrate 
how  to  apply  these  to  the  present  model.  We  found  through  numerical  experiments  that  the 
WENO  method  is  vastly  more  efficient  than  the  low  order  Lax-Friedrichs  scheme,  yet  both 
methods  converge  to  the  same  solution  of  the  physical  model.  It  is  especially  interesting  to 
observe  the  small  staircases  in  the  solution  which  are  completely  missed  out,  because  of  the 
numerical  viscosity,  if  a  lower  order  method  is  used  without  a  sufficiently  refined  mesh.  To 
demonstrate  the  applicability  of  this  new,  efficient  numerical  tool,  we  study  the  multi-class 
model  under  different  parameter  regimes  and  traffic  stream  models.  We  consider  also  the 
convergence  of  the  multi-class  LWR  model  when  the  number  of  classes  goes  to  infinity.  We 
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1.  INTRODUCTION 

Lighthill  and  Whitham  (1955)  and  Richards  (1956)  independently  proposed  a  simple 
continuum  model,  now  known  as  the  LWR  model,  to  describe  the  characteristics  of  traffic 
flow.  In  this  model,  a  traffic  stream  model  (a  relationship  between  the  traffic  state  variables 
of  flow,  speed  and  density,  e.g.  Greenshields  (1934))  is  supplemented  by  the  continuity 
equation  of  vehicles,  and  the  resultant  partial  differential  equation  presumably  could  be 
solved  to  obtain  the  density  as  a  function  of  space  and  time.  For  a  specific  form  of 
Greenshields’  traffic  stream  model,  the  solution  can  be  obtained  analytically  (Wong  and 
Wong,  2002).  Although  aiming  to  provide  a  coarse  representation  of  traffic  behavior,  the 
LWR  model  is  capable  of  reproducing  qualitatively  a  remarkable  amount  of  real  traffic 
phenomena  such  as  shock  formation.  Nevertheless,  there  are  still  some  puzzling  traffic 
phenomena  observed  on  the  highway,  such  as  the  two-capacity  or  reverse-lambda  state  in  the 
fundamental  diagram,  hysteresis  of  traffic  flow  and  platoon  dispersion,  that  this  simple  LWR 
model  cannot  address  or  explain. 

Recently,  a  multi-class  model,  extended  from  the  LWR  model,  with  heterogeneous  drivers 
was  formulated  (MCLWR  model)  (Wong  and  Wong,  2001).  Let  there  be  M  classes  of  road 
users  with  different  speed  choice  behaviors  in  response  to  the  same  traffic  density  when 
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traveling  on  a  highway  section.  It  means  that  for  a  given  total  density,  there  exists  a 
distribution  of  equilibrium  speeds  by  different  user  classes.  It  is  expected  that  the  variation 
around  the  mean  speed  (averaged  over  all  user  classes)  decreases  when  traffic  density 
increases,  due  to  the  tighter  interactions  between  road  users.  Let  qm(x,t),  km(x,t)  and 
um(x,t )  be,  respectively,  the  flow,  density  and  speed  of  user  class  m  in  the  space-time 
domain.  The  total  density  on  a  highway  section  can  then  be  obtained  as 

M 

k(x,t)=YJkm{x,t).  (1) 

m= 1 


The  flow,  density  and  speed  variables  of  a  particular  class  are  subject  to  the  following 
definitional  relationship, 


qm  (x,  t )  =  um  (x,  t )  •  km  (x,  t),  Vm  =  1,2, . . . ,  M  . 


(2) 


From  the  law  of  conservation  of  vehicles,  each  user  class  should  satisfy  the  following 
continuity  equation, 


dt 


dx 


(3) 


which  describes  the  conservation  of  vehicles  at  any  location  at  any  time  along  a 
topographically  homogeneous  highway  section  without  intermediate  entrances  or  exits. 
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The  core  of  the  present  extension  is  to  assume  that  the  choice  of  speed  of  a  particular  user 
class  is  not  only  affected  by  the  presence  of  this  user  class,  but  also  by  all  other  user  classes 
on  the  highway.  A  general  fonn  of  speed-density  relationship  can  be  written  as 


um  (x, t )  =  Um  (kx,k2,...,kM),  Mm  =  1,2, ...,M  . 


(4) 


For  the  isotropic  case,  the  above  relationship  would  take  a  simpler  functional  form  as 

um (x, t )  =  Um (k),  Mm  =  1,2, ...,M  ,  (5) 


where  k  is  the  total  density  determined  by  equation  (1). 

Combining  the  above  equations,  the  problem  can  be  fonnulated  into  a  set  of  partial 
differential  equations, 


where 


dk 


!  yc 

*  h 


dkn(x,t ) 

mn\ u’ 
OX 


Mm  =  1,2,  ...,M 


=  U„8„„+k„  dUm 


ok, 


Mm,n  =  1,2, ...,M , 


(6) 


(7) 


is  the  kinematic  wave  speed  of  user  class  m  in  response  to  the  presence  of  class  n  users,  and 
8mn  =1  if  m  =  n;  and  8mn  =  0  if  m  ^  n.  Note  that  the  problem  stipulated  in  equation  (6) 

reduces  to  the  original  LWR  model  when  M  =  1  (i.e.  homogeneous  users).  The  problem 
becomes  one  of  solving  the  set  of  differential  equations  (6),  or  better  still  the  conservation 
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form  (3)  with  qm  defined  by  equation  (2)  and  um  defined  by  equation  (5),  subject  to  certain 


initial  spatial  and  time  boundary  conditions.  Although  the  problem  formulation  is  straight 
forward,  it  was  found  that  the  model  is  capable  of  producing  the  desired  properties  of  a 
macroscopic  traffic  flow  model  and  it  explains  many  puzzling  phenomena,  such  as  the  two- 
capacity  or  reverse-lambda  state,  hysteresis,  and  platoon  dispersion,  but  it  would  not  be 
subject  to  other  deficiencies  such  as  wrong-way  travel  (Daganzo,  1995). 

In  Wong  and  Wong  (2001),  the  MCLWR  model  was  solved  by  a  first  order  Lax-Friedrichs 
finite  difference  scheme.  Although  this  finite  difference  scheme  is  commonly  used  to  solve 
the  original  LWR  model  (Lebacque,  1984;  Michalopoulos  et  ah,  1984),  it  is  argued  that  this 
first  order  Lax-Friedrichs  scheme  may  produce  smeared  solutions  near  discontinuities  due  to 
excessive  numerical  viscosity.  The  effect  of  numerical  viscosity  will  diminish  with  mesh 
refinement,  but  it  will  be  very  costly  to  solve  a  very  refined  mesh.  More  recently,  Lebacque 
(1996)  successfully  applied  the  Godunov  scheme,  introduced  by  Godunov  (1959),  to  solve 
the  LWR  model.  The  Godunov  scheme  is  subject  to  smaller  numerical  viscosity,  but  it 
requires  a  Riemann  solver  as  its  building  block,  which  is  very  difficult,  if  not  impossible,  to 
develop  for  the  MCLWR  model.  This  is  because  the  multi-class  model  does  not  seem  to  be 
either  genuinely  nonlinear  or  linearly  degenerate  (LeVeque,  1992).  Nevertheless,  it  is 
important  to  note  that,  even  though  for  first  order  methods  the  Godunov  scheme  is  more 
accurate  than  the  Lax-Friedrichs  scheme,  this  difference  diminishes  dramatically  when  higher 
order  schemes  are  considered  (Shu,  1998).  Both  Godunov  and  Lax-Friedrichs  schemes 
converge  to  the  same  physical  solution  of  the  model  with  a  sufficiently  refined  mesh.  This 
can  be  proved  for  the  scalar  and  some  system  cases,  and  can  be  observed  for  more  complex 
systems  (LeVeque,  1992;  Shu,  1998). 
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This  paper  presents  the  solution  of  the  MCLWR  model  by  a  weighted  essentially  non- 
oscillatory  (WENO)  scheme  (Jiang  and  Shu,  1996).  The  WENO  scheme  is  a  very  robust 
numerical  scheme  and  is  found  to  be  very  useful  in  computational  fluid  dynamics  as  well  as 
in  other  applications.  The  numerical  results  from  WENO  are  compared  with  those  obtained 
from  the  first  order  Lax -Friedrichs  method.  In  the  special  case  when  all  the  eigenvalues  of  the 
kinematic  wave  matrix  of  the  system  (6)  are  positive,  the  Godunov  solver  becomes  the  simple 
upwind  solver.  In  this  special  case  we  have  verified  that  the  first  order  Godunov  solver 
converges  faster  than  the  first  order  Lax-Friedrichs  solver,  but  slower  than  the  fifth  order 
WENO  solver,  while  all  three  converge  to  the  same  physical  solution.  In  Section  2,  some 
analytical  properties  of  the  MCLWR  model  are  presented  for  a  2-class  model.  The  WENO 
scheme  for  the  MCLWR  is  given  in  Section  3.  Section  4  compares  the  convergence 
characteristics  of  the  numerical  schemes,  shows  the  numerical  solutions  for  different 
congestion  regimes,  and  studies  the  asymptotic  case  of  an  infinite  number  of  classes. 

2.  SOME  ANALYTICAL  PROPERTIES  OF  THE  MCLWR  MODEL 

2,1  Hyperbolicity  of  the  system 

Let  k  =  Col (kl,k2,-..,kM)  be  a  column  vector  containing  all  M  density  variables  in  the 
MCLWR  model.  For  a  smooth  solution  k  (meaning  that  k  has  at  least  first  order  continuous 
derivatives  in  x  and  t),  the  set  of  partial  differential  equations  (PDEs)  (3)  can  be  rewritten  as 


5k 

dt 


.  n  .  5k 
+  A(k)— 

ox 


=  0, 


(8) 


6 


where  A(k)  =  Vkq(k)  is  a  kinematic  wave  matrix  containing  all  the  elements  of  cm„  in 


equation  (7),  and  q  =  Col  (q\,q2,---,qM)  is  a  column  vector  of  the  flow  fluxes  of  all  classes. 
However,  when  the  solution  k  becomes  discontinuous  (containing  shocks  or  other 
discontinuities),  the  two  systems  (3)  and  (8)  are  not  equivalent.  Thus  to  be  on  the  safe  side, 
one  should  always  use  conservative  schemes  to  solve  (3)  directly. 

The  system  (3)  is  called  hyperbolic  if  the  eigenvalues  of  the  kinematic  wave  matrix  A(k)  are 
all  real  and  there  is  a  complete,  linearly  independent  set  of  eigenvectors.  Hyperbolic  systems 
are  mathematically  well-posed,  meaning  that  their  solutions  depend  continuously  on  the 
initial  conditions.  This  can  be  proved  for  the  linear  case  and  also  for  some  nonlinear  cases. 
An  important  issue  in  verifying  the  reasonableness  of  a  model  is  to  check  if  it  is  hyperbolic. 
We  confirm  that  the  system  (3)  is  hyperbolic  for  the  practical  choices  of  traffic  stream  models 
and  their  parameters.  This  verification  is  perfonned  analytically  for  the  2-class  case  and 
numerically  for  the  general  M-class  case  to  be  discussed  in  Section  4. 

For  the  2-class  case  (M=  2),  the  kinematic  wave  matrix  is  a  2  x  2  matrix  given  by 


A  = 


Ux(k)  +  kxU[(k) 
k2U'2{k) 


k\U\(k) 

U2(k)  +  k2U'2(k)  J’ 


(9) 


where  U\(k)  and  U2(k)  are  defined  in  equation  (5).  The  two  eigenvalues  of  the  kinematic 
wave  matrix  are  thus  given  by 

>-1,2  =(ul(k)  +  klU[(k)  +  U2(k)  +  k2U'2(k)±y[D)/2,  (10) 

where 
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D  =  ((Ul  (k)  +  kxU[  (k)) -  (U 2  (k)  +  k2U'2  (k))f  +  4k} k2U ,  (Ar)C/2  (*)  • 


(ID 


Clearly  D  >  0  which  implies  that  the  two  eigenvalues  are  both  real  and  distinct.  Thus  the 
system  is  hyperbolic. 

We  could  obtain  analytical  fonnulas  for  the  eigenvalues  of  the  kinematic  wave  matrix  for  the 
3-class  or  even  the  4-class  case  (M  =  3  or  4),  however  the  formulas  are  quite  complex  and  it 
is  not  easy  to  see  whether  the  eigenvalues  are  always  real.  At  any  rate,  this  approach  would 
not  work  for  the  multi-class  case  with  M  >  4,  as  no  analytical  formulas  for  the  eigenvalues 
would  be  available. 

We  thus  resort  to  implementing  a  numerical  eigenvalue  solver  to  verify  a  posteriori  that  the 
eigenvalues  are  always  real  for  all  the  test  cases  in  Wong  and  Wong  (2001)  and  in  this  paper. 
Indeed  through  all  the  numerical  tests,  non-real  eigenvalues  never  appear  for  the  kinematic 
wave  matrix.  Although  this  is  not  a  rigorous  proof  that  the  MCLWR  model  is  always 
hyperbolic,  it  at  least  gives  validity  to  the  numerical  experiments  in  Wong  and  Wong  (2001) 
and  in  this  paper  since  the  models  under  all  these  cases  are  hyperbolic. 

2.2  First  order  traveling  waves 

In  this  section,  we  apply  a  linearization  approach  to  demonstrate  the  traveling  wave 
properties  of  the  MCLWR  model  (Whitham,  1974).  To  simplify  the  analysis,  we  also 
consider  the  simple  2-class  system  and  assume  a  modified  Greenshields’  fonn  of  traffic 
stream  model, 
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(12) 


s 

II 

s 

1  k\+k2 

Jr . 

and  u2  =  u  y2 

1  kx+k2 

If . 

K  Jam  ) 

\  Jam  y 

where,  for  Class  1  and  Class  2  traffic,  u\  and  112  are  the  traffic  speeds,  k\  and  C  are  the 
densities,  and  uj \  and  u/2  are  the  free-flowing  speeds,  respectively,  while  k]Am  is  the  jam 
density  of  the  highway.  The  system  of  differential  equations  can  be  written  as 

dkx  d(uxkx )  ,  dk2  d{u2k2 

dt  dx  dt  dx 

For  small  perturbations  of  densities,  r  and  w,  around  the  steady  state  densities  kx  and  k2 
respectively,  we  can  write 


=  0. 


(13) 


kx  =  kx  +  r  and  k2  =k2+w .  (14) 

Substituting  equation  (14)  into  equation  (13)  and  neglecting  higher  order  terms,  we  can  show 
that 


dr 

—  + 

dt 


Vi 


k\ 


( k 


jam 


jam 


—  —  dr  Ufxkx  dw 

2k\-k2)—  =  -LL-L— 
dx  A^jam 


and 


dt  k- 

^jam 


jam 


r  T  dw  uf2  T  dr 
■  2A't  -  kx )  —  =  — —  k2  — . 

dx  Arjam  “  dx 


(15) 


(16) 


Eliminating  w  from  equations  (15,  16),  we  have 
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(17) 


where 


d  <9..  3  d  .  k\k2u  nu  f2  d2  r 

(—  +  ®  1  — )(—  +  ®2  —)r  = - — - — , 

dt  ox  dt  ox  kr  dx~ 

^jam 


k]U  ry  knU  f-) 

coj  =iq - : —  and  (d2  —u2  — 


ljam 


k  i 


jam 


(18) 


In  equation  (18),  the  left-hand  side  is  a  wave  operator,  which  indicates  that  there  are  two  first 
order  traveling  waves  of  speeds  coj  and  co  2  in  the  traffic  stream.  It  is  interesting  and  important 
to  note  that  these  class-characterized  waves  always  travel  more  slowly  than  the  fastest  vehicle 
in  the  traffic  stream.  This  linearization  approach  can  also  be  generalized  to  any  number  of 
classes.  Also  note  that  when  substituting  the  modified  Greenshields’  form  of  traffic  stream 
model  (12)  into  equation  (10),  the  eigenvalues  are  identical  to  the  traveling  wave  speeds 
shown  in  equation  (18)  for  this  2 -class  case,  when  one  of  the  steady  state  densities  kx  or  k2 
is  equal  to  zero. 


3.  WENO  NUMERICAL  SCHEME 


In  recent  years  many  high  order,  high  resolution,  numerical  methods  have  been  developed  in 
the  literature  to  solve  a  system  of  partial  differential  equations  (PDEs).  The  main  applications 
are  in  computational  fluid  dynamics,  but  there  are  also  applications  in  other  physical  and 
engineering  areas.  In  this  paper  we  apply  the  high  order  finite  difference  scheme,  weighted 
essentially  non-oscillatory  (WENO)  scheme  (Jiang  and  Shu  1996;  Balsara  and  Shu,  2000),  to 
solve  the  MCLWR  system  (3).  In  particular,  the  fifth  order  WENO  scheme  in  Jiang  and  Shu 
(1996)  is  used.  The  numerical  procedure  is  summarized  in  this  section.  These  numerical 
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methods  are  found  to  be  very  useful  because  of  their  simultaneous  high  order  accuracy  and 
non-oscillatory  property  in  the  presence  of  shocks  and  other  discontinuities  or  sharp  gradient 
regions  in  the  solution,  or,  in  general,  for  convection  dominated  problems.  For  more  details  of 
such  methods,  see  Shu  (1998)  and  Shu  (2002). 

We  now  describe  the  computational  procedure  of  the  fifth  order  WENO  scheme.  Spatial 
discretization  is  discussed  first.  We  start  from  the  simple  case  of  a  scalar  equation  (3),  i.e.  a 

dq(k) 

l-class  model,  and  assume  ^  >  0 ,  i.e.  the  wind  direction  is  positive.  More  general  cases 

will  be  described  later.  The  computational  domain  is  discretized  into  a  unifonn  mesh  of  J  grid 
points: 


Xj  =  j'Ax  ;  j  =  1,2,...,  J ,  (19) 

where  Ax  is  the  uniform  mesh  size  on  the  spatial  axis.  A  smooth  non-uniform  mesh  could 
also  be  used  to  concentrate  grid  points  near  certain  regions  to  obtain  better  resolution.  A 
conservative  numerical  approximation  kj(t)  to  the  exact  solution  k(xj,t )  of  (3)  satisfies  the 
following  ordinary  differential  equation  (ODE)  system: 

d  k  (t)  i  /  \ 

-^r  +  -fc+l/2-^l/2)=0,  (20) 

where  cjj+i/2  is  called  the  numerical  flux,  whose  design  is  the  key  ingredient  for  a  successful 
scheme.  For  the  fifth  order  WENO  scheme,  the  numerical  flux  qJ+ 1/2  is  defined  as  follows: 
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(21) 


r,  _n=(!)  ,  n  ;=(2)  ,n=(3) 

<7 7+1/2  -  ult7 7+1/2  +  U2<7 7+1/2  +  U3<7 y+1/2  ’ 


where  q^fXll  are  three  third  order  fluxes  on  three  different  stencils  given  by 


-0)  1  7  11 

^7+1/2  =3^7-2  -pj-l  +J4j’ 

~(2)  1  5  1 

tf)+i/2  =_6^-t  +  g' qj  +  3^+t’ 

-(3)  1  ,  5  1 

?}+i/2  =  3^-  +  6^+1-  6^-+2’ 


and  the  nonlinear  weights  u/(  are  given  by 


u 


p 


Y  / 

(e  +  P/)2’ 


with  the  linear  weights  y/  given  by 


Yi  = 


(22a) 

(22b) 

(22c) 


(23) 


(24) 


and  the  smoothness  indicators  (3/  given  by 


Pi  ={2^,-2  -'iqj-x  +  qj)2  +\(qj-2  -4qj-i  +  ^j)2 ■ 


(25) 
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(26) 
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P2  =  —(qj-i-2qJ+qJ+ir  +-(<//-,  -<7y+i)  , 


Ps  =  ]|<A/  - 2q j+\  +qj+2 )2  +  ^(3^  -4^y+l  +^+2)2  •  (27) 


Finally,  s  is  a  parameter  to  prevent  the  denominator  from  becoming  0  and  is  fixed  at 

s  =  10  6  in  all  the  computations  in  this  paper.  The  choice  of  s  does  not  affect  accuracy:  the 
numerical  errors  can  go  much  lower  than  s ,  reaching  machine  zeros  (around  10'  for  double 
precision).  Notice  that  we  have  used  the  short  notation  qj  to  denote  q(kj(t)).  We  remark  that 
the  stencil  for  the  scheme  is  biased  to  the  left  because  of  the  positive  wind  direction. 

This  finishes  the  description  of  the  fifth  order  finite  difference  WENO  scheme  (Jiang  and  Shu 
1996)  for  the  scalar  equation  with  a  positive  wind  direction.  As  we  can  see,  the  algorithm  is 
actually  quite  straight  forward  and  there  are  no  parameters  to  be  tuned  in  the  scheme.  The 
main  reason  that  it  works  well,  both  for  smooth  solutions  and  for  solutions  containing  shocks 
or  other  discontinuities  or  high  gradient  regions,  is  that  the  nonlinear  weights,  determined  by 
the  smoothness  indicators,  are  automatically  adjusting  themselves,  based  on  the  numerical 
solution,  to  use  the  locally  smoothest  information  given  by  the  solution.  We  refer  to  (Jiang 
and  Shu  1996;  Balsara  and  Shu  2000;  Shu  1998;  and  Shu,  2002)  for  accuracy  tests  and 
computational  fluid  dynamics  simulations  using  this  method. 

If  the  “wind  direction”  dq^  -  <  0 ,  the  procedure  for  computing  the  numerical  flux  qj+ 1/2  is  a 
mirror  image  with  respect  to  the  point  xj+  m,  of  what  is  described  above.  The  stencil  would 

da(k) 

then  be  biased  to  the  right.  If  — changes  sign,  we  will  use  smooth  flux  splitting 
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q(k)  =  q+(k)  +  q  ( k ) 


(28) 


where  dq ^  >  0  and  dq ^  <  0 ,  and  apply  the  above  procedures  separately  on  each  one  of 

them.  There  are  many  choices  of  such  flux  splitting,  however  the  most  popular  one  is  the 
Lax-Friedrichs  flux  splitting  where 


q±(k)  =  ^(q(k)±ak). 


(29) 


with  a  =  max  k 


Sq(k) 

8k 


For  hyperbolic  systems  of  conservation  laws  (3),  the  eigenvalues  of  A(k)  are  all  real: 

kx{k)<..<XM{k).  (30) 

A  safe  but  rather  expensive  way  to  generalize  scalar  schemes  to  such  system  cases  is  to 
utilize  local  characteristic  decompositions,  see  Shu  (1998)  for  details.  However,  such  a 
procedure  depends  on  the  explicit  fonnulas  for  the  eigenvectors  of  A(k) ,  which  are  not  easy 
to  obtain  for  the  MCLWR  model  when  M  >  3  .  We  have  thus  adopted  a  simpler,  component¬ 
wise  generalization,  namely  using  the  Lax-Friedrichs  flux  splitting  (29)  for  each  equation  in 
the  system,  with  a  common  a .  Ideally,  a  should  be  chosen  as  the  largest  (absolute  value) 
eigenvalue  in  (30),  however  closed  fonn  formulas  for  the  eigenvalues  (30)  are  also  difficult 
to  obtain  for  the  MCLWR  model  when  M  >  3  .  We  have  thus  chosen  a  as 
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a  =  ma x(iq um  ), 


(31) 


where  the  speeds  um  are  defined  in  equation  (5).  Numerical  experiments  in  the  next  section 
indicate  that  this  works  well  for  the  MCLWR  model. 

For  the  time  discretization,  the  computational  domain  is  discretized  into  a  mesh  of  N  grid 
points:  tn  =  t"~l  +  At'1;  n  =  1,2,..., TV  ,  where  At"  is  the  mesh  size  on  the  time  axis.  We  adopt 
the  third  order  TVD  Runge-Kutta  method  (Shu  and  Osher  1988): 

kil)  =kn  +AtnL(kn,tn), 

k(2)  =-k"  +-k{1)  +-At"L(km,t"  +Atn), 

4  4  4 

k (3)  =-k"  +-k{2)  +-A  t"L(k{2\t"  +-A  tn), 

3  3  3  2 

where  L  is  the  approximation  to  the  spatial  derivatives: 

(33) 

OX 

established  by  the  WENO  procedure  outlined  above.  This  time  discretization  is  proved  stable 
if  the  first  order  Euler  forward  time  stepping  of  the  spatial  operator  is  stable  (Shu  and  Osher, 
1988;  Gottlieb  et  al.,  2001).  Note  that  this  time  discretization  is  very  simple  and  consists  of 
convex  combinations  of  three  first  order  Euler  forward  steps.  A  CFL  condition  is  needed  for 
stability: 


(32a) 

(32b) 

(32c) 
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(34) 


a 


A  tn 

n - <  CFL  . 

Ax 


where  a"  should  be  taken  as  the  largest  (absolute  value)  eigenvalue  in  (30)  for  time  level  n, 
but  for  the  MCLWR  model  it  is  taken  as  that  in  (31)  instead.  CFL  should  be  less  than  one  for 
stability  and  in  our  computation  it  is  taken  as  0.6. 

4.  COMPUTATIONAL  EXPERIMENTS 

In  this  section  we  present  our  computational  experiments  of  the  MCLWR  model  using  the 
high  order  WENO  scheme,  and  compare  the  results  with  other  numerical  schemes. 

4,1  Traveling  wave  speeds 

We  start  with  a  2-class  model  with  the  modified  Greenshields’  fonn  of  traffic  stream  model 
as  in  equation  (12).  Consider  a  highway  2  km  long  with  an  initial  platoon  of  maximum 
density  40  veh/km  as  shown  in  Figure  1 .  The  left  boundary  has  no  inflow  (density  equals  0) 
for  all  time,  and  the  right  boundary  is  a  free  outflow  (Neumann  boundary  condition).  The 
free-flowing  speeds  of  Class  1  and  Class  2  drivers  are  60  km/h  and  120  km/h  respectively. 
We  assume  an  equal  distribution  of  drivers  in  the  platoon.  The  jam  density  of  the  highway  is 
200  veh/km.  The  dispersion  of  the  platoon  at  time  0.01  hours  is  shown  in  Figure  2.  It  is 
interesting  to  note  that  the  solution  forms  two  uniform  density  platforms  or  staircases,  (to  be 
discussed  later),  each  of  which  contains  only  a  single  class  of  driver.  The  widths  of  the 
platforms  are  marked  by  points  A,  B,  C  and  D  in  the  figure,  with  the  platfonn  A-B  containing 
Class  1  drivers  only,  and  the  platform  C-D  containing  Class  2  drivers  only.  In  Figure  3,  we 
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show  the  speed  trajectories  of  these  points,  A,  B,  C  and  D,  which  numerically  measure  the 
wave  speeds,  together  with  the  values  obtained  from  the  eigenvalue  fonnulas  (10)  and  those 
given  by  the  linearized  formulas  (18).  The  result  shows  that  the  linearized  wave  speeds  give 
reasonable  predictions  of  actual  nonlinear  wave  speeds  in  this  case. 

4.2  Convergence  study  of  the  numerical  methods 

For  a  nonlinear  system  (3)  with  possible  shocks  and  other  discontinuities,  it  is  not  possible  to 
prove  mathematically  that  the  WENO  scheme,  or  any  other  scheme,  converges.  However, 
experience  from  computational  fluid  dynamics  indicates  that  the  WENO  scheme  is  very 
robust  and  always  converges  for  hyperbolic  systems.  We  would  like  to  verify  through 
numerical  experiments  the  convergence  of  the  WENO  scheme  for  the  MCLWR  model  in  this 
subsection. 

For  this  purpose  we  take  Experiment  2  of  Wong  and  Wong  (2001)  as  our  test  case.  Other  test 
cases  have  also  been  experimented  with,  yielding  similar  results.  In  this  experiment,  we 
consider  the  same  highway  section  and  initial  density  platoon  as  that  in  Section  4.1,  but  with 
the  number  of  driver  classes  increased  to  M  =  9.  The  traffic  stream  model  takes  the  modified 
Drake’s  form  (Drake  et  ah,  1967)  as 

=Um(k)  =  ufmexp(-(k/k0)2 /2\  m  =1,2,...,M  .  (35) 

The  free  flowing  speeds  Ufm  of  these  drivers  are  taken  as  60.0,  67.5,  75.0,...,  120.0  km/h,  and 
the  optimal  density  ko  =  50  veh/km.  The  distribution  in  density  for  these  user  classes  is  given 
by  Figure  4.  The  left  boundary  has  no  inflow  (density  equals  0)  for  all  time,  and  the  right 
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boundary  is  a  free  outflow  (Neumann  boundary  condition).  We  plot  the  density  at  t  =  0.015 
hours  to  verify  numerical  convergence. 

We  first  plot,  in  Figure  5,  top  left,  the  total  density  computed  with  the  first  order  Lax- 
Friedrichs  scheme  (described  in  detail  in  Wong  and  Wong  (2001))  using  6400  grid  points 
(solid  line),  versus  that  computed  with  the  fifth  order  WENO  scheme  using  100  grid  points 
(circles).  They  overlay  each  other  quite  well,  indicating  two  things: 

1.  The  resolution  of  the  first  order  Lax -Friedrichs  scheme  with  6400  points  is  similar  to  that 
of  the  fifth  order  WENO  scheme  with  100  points.  Thus  the  high  order  WENO  scheme  is 
vastly  more  efficient  than  the  first  order  Lax-Friedrichs  scheme  for  this  test  case.  This 
conclusion  is  also  valid  in  general. 

2.  If  one  uses  only  a  first  order  scheme  to  compute,  one  might  decide  prematurely  that  this  is 
a  convergent  solution,  since  6400  points  make  a  very  refined  mesh. 

In  fact,  the  solution  in  Figure  5,  top  left,  is  not  a  convergent  one  numerically,  although  the 
solution  is  good  enough  to  demonstrate  the  physical  characteristics  of  the  traffic  model. 
Nevertheless,  in  this  paper,  we  study  the  numerical  convergence  characteristics  of  the  traffic 
model  in  greater  detail.  In  Figure  5,  top  right,  we  plot  the  WENO  solutions  using  100  points 
(dash-dot  line),  400  points  (dashed  line)  and  1600  points  (solid  line).  We  can  see  that  the 
solution  has  observable  differences  for  all  these  grids  (the  difference  between  the  400  points 
and  1600  points  results  is  small  but  still  noticeable,  especially  when  enlarged  near  the 
staircases).  In  particular,  notice  the  small  staircases  in  the  increasing  part  of  the  solution. 
There  are  9  such  small  staircases,  clearly  related  to  the  9  user  classes.  The  coarse  mesh  (100 
points)  WENO  solution  and  most  of  the  first  order  Lax-Friedrichs  solutions  completely  miss 
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these  staircases  because  of  the  excessive  numerical  dissipation.  To  verify  that  the  numerical 
solution  is  indeed  convergent,  we  have  also  computed  it  using  WENO  with  3200  grid  points. 
The  solution  (not  shown)  completely  overlays  that  with  1600  grid  points,  indicating  that  the 
WENO  solution  with  1600  grid  points  can  be  considered  a  numerically  convergent  solution. 

In  Figure  5,  bottom  left,  we  plot  the  first  order  Lax-Friedrichs  solutions  using  400  points 
(dash-double  dot  line),  1600  points  (dash-dot  line),  6400  points  (dashed  line)  and  25600 
points  (solid  line).  We  can  see  that  the  solution  does  eventually  converge  with  grid 
refinements,  however  such  convergence  is  very  slow  and  one  needs  a  huge  number  of  grid 
points  (in  this  case  25600  points)  to  see  the  staircases.  To  convince  the  reader  that  both  the 
WENO  scheme  and  the  Lax-Friedrichs  scheme  converge  to  the  same  solution,  in  Figure  5, 
bottom  right,  we  plot  the  WENO  solution  using  400  points  (circles)  and  the  Lax-Friedrichs 
solution  using  25600  points  (solid  line).  They  overlay  each  other  quite  well,  both  showing  the 
small  staircases. 

In  Figure  6,  left,  we  plot  the  total  density  as  a  function  of  spatial  location,  for  various  times,  t 
=  0,  0.005,  0.010,  ...,  0.025  hours.  We  can  clearly  see  the  evolution  of  the  dispersion  of  the 
back  of  the  platoon  and  the  appearance  of  staircases.  In  Figure  6,  right,  we  plot  the  flow 
(defined  as  the  sum  of  the  fluxes  in  (2)  over  all  M  classes)  as  a  function  of  time,  t,  at  various 
spatial  locations,  x  =  0.2,  0.4,  0.6,...,  2.0  km.  We  can  clearly  see  that  the  small  staircases  are 
also  present  in  these  flow  plots.  The  results  in  Figure  6  are  obtained  using  WENO  with  1600 
grid  points,  which  overlays  well  the  WENO  results  using  3200  grid  points  (not  shown  here), 
indicating  that  they  are  reliable,  numerically  convergent  solutions  of  the  9-class  model. 
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To  assure  the  reader  that  these  staircases  are  not  numerical  artifacts  of  the  WENO  schemes, 
we  plot  in  Figure  7  the  total  density  as  a  function  of  spatial  location,  for  various  times,  t  =  0, 
0.005,  0.010, ...,  0.030  hours,  using  the  9-class  model  but  with  Ufm  =  90  km/h  for  all  9  classes. 
Clearly  this  returns  to  a  single  user  model  and  the  solution  is  now  free  from  the  staircases. 
The  results  in  Figure  7  are  obtained  using  WENO  schemes  with  1600  grid  points,  which 
overlay  the  WENO  results  using  3200  grid  points  (not  shown  here),  indicating  that  they  are 
reliable,  numerically  convergent  solutions.  Note  that  a  single  user  model  is  the  original  LWR 
model,  which  shows  no  dispersion  behavior,  as  revealed  by  the  figure. 

In  summary,  in  this  subsection  we  have  shown  that 

1.  Using  the  fifth  order  WENO  scheme  is  vastly  more  efficient  than  using  the  first  order 
Lax-Friedrichs  scheme,  saving,  by  a  factor  of  64,  on  the  number  of  mesh  points  needed  to 
reach  the  same  resolution; 

2.  One  must  be  very  careful  in  performing  the  grid  refinement  study  to  verify  numerical 
convergence,  for  otherwise  one  might  miss  some  very  important  solution  features,  such  as 
the  staircases,  which  might  otherwise  be  completely  obscured  by  numerical  dissipation; 

3.  Both  the  high  order  WENO  scheme  and  the  low  order  Lax-Friedrichs  scheme  eventually 
converge  to  the  same  physical  solution  with  grid  refinements. 

Finally  in  this  subsection,  we  point  out  that  we  have  verified  a  posteriori,  by  a  numerical 
eigenvalue  solver,  that  all  the  eigenvalues  of  the  kinematic  wave  matrix  stay  non-negative 
during  the  time  evolution.  We  remark  that  this  is  true  only  for  this  case  and  not  in  general  for 
cases  in  the  next  subsection.  For  this  special  case,  it  is  straight  forward  to  write  out  the  first 
order  Godunov  scheme,  which  coincides  with  the  simple  upwind  scheme  (using  backward 
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difference  to  approximate  the  spatial  derivatives).  In  Figure  8,  we  show  the  convergence 
history  as  well  as  the  density  and  flow  graphs  computed  by  the  first  order  Godunov  scheme. 
We  can  clearly  see  that  the  resolution  of  the  first  order  Godunov  scheme  is  better  than  that  of 
the  first  order  Lax-Friedrichs  scheme  but  worse  than  that  of  the  fifth  order  WENO  scheme. 
Moreover,  the  first  order  Godunov  scheme  converges  to  the  same  solution  as  the  other  two 
schemes,  when  we  overlay  the  solutions  (not  shown  here). 

The  first  order  Godunov  scheme  has  merit  in  that  it  is  computationally  fast.  However,  the 
clocked  times  (shown  in  Table  1)  for  both  the  fifth  order  WENO  code  and  the  first  order 
upwind  code  (Godunov  in  the  special  case  of  all  positive  eigenvalues),  which  achieve  the 
same  resolution,  still  favor  the  WENO  scheme.  Clearly  the  advantage  of  the  simple  and  fast 
computation  of  the  first  order  Godunov  scheme  is  offset  by  the  use  of  more  grid  points  to 
achieve  higher  accuracy.  We  also  point  out  again  that  the  Godunov  scheme  for  this  system  is 
very  difficult,  if  not  impossible,  to  obtain  when  the  eigenvalues  change  sign. 

4.3  Numerical  experiments  for  different  congestion  regimes  and  traffic  stream 
models 

In  this  subsection,  we  perform  more  numerical  experiments,  using  the  WENO  scheme  with 
1600  grid  points  (which  gives  numerically  convergent  solutions),  with  different  model 
characteristics  for  the  MCLWR  model.  We  still  use  the  9-class  model  with  the  density 
distribution  given  by  Figure  4,  but  we  consider  the  following  four  cases: 
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1 .  The  initial  density  distribution  represents  a  platoon  in  the  non-congested  regime,  as  given 
in  Figure  1.  The  traffic  stream  model  takes  the  modified  Drake's  form  (35).  This  case  has 
already  been  considered  in  the  previous  subsection; 

2.  The  initial  density  distribution  represents  a  platoon  in  the  congested  regime,  as  given  in 
Figure  9,  which  has  a  higher  maximum  initial  density  value  of  120  veh/km.  The  traffic 
stream  model  is  identical  to  that  in  Case  1 ; 

3.  The  initial  density  distribution  represents  a  platoon  in  the  non-congested  regime,  as  given 
in  Figure  1.  The  traffic  stream  model  takes  the  modified  Greenshields’  fonn  (36)  with  a 
jam  density  of  Ajam  =  200  veh/km, 

~  U m  (^-)  —  M  fm  (l —  k  /  Ajam  j,  ft l  —  1,2,  —  ,M  ,  (36) 

4.  The  initial  density  distribution  represents  a  platoon  in  the  congested  regime,  as  given  in 
Figure  9.  The  traffic  stream  model  is  identical  to  that  in  Case  3. 

The  density  versus  distance  plots  for  various  times,  and  the  flow  versus  time  plots  for  various 
spatial  locations,  are  given  in  Figures  6,  10,  1 1  and  12,  respectively,  for  these  cases.  It  is  clear 
from  these  figures  that  dispersion  at  the  tail  of  the  platoon  is  limited  until  its  density  value  has 
dropped  to  near  or  below  certain  critical  values.  These  values  are  the  optimal  densities  given 
by  the  Drake’s  and  Greenshields’  traffic  stream  models  used.  They  are  related  to  the  capacity 
or  maximum  flow  of  the  highway  being  analyzed.  A  given  highway  is  described  as  operating 
at  a  congested  or  non-congested  state  when  the  traffic  stream  is  flowing  above  or  below  the 
capacity  of  the  highway,  respectively. 
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We  first  look  at  Figure  6  for  the  modified  Drake’s  model.  The  optimal  density,  Ay  is  set  to  be 
50  veh/km  for  this  model.  Initially  the  density  of  the  platoon  is  below  the  optimal  density  and 
therefore  the  highway  is  not  congested.  As  a  result  vehicles  in  the  9  different  classes  are  free 
to  overtake  and  move  as  desired,  and  the  platoon  disperses  as  depicted  in  the  figure.  The 
second  case  for  Drake’s  model  is  an  initial  platoon  of  density  120  veh/km,  much  greater  than 
the  highway’s  optimal  density.  This  time  the  highway  is  operating  in  a  congested  state  and 
this  results  in  a  non-dispersed  tail  of  the  platoon  because  overtaking  is  limited  in  such  a 
congested  state.  The  front  of  the  platoon  still  can  disperse  because  the  downstream  end  is 
empty.  Vehicles  at  the  tail  of  the  platoon  however  have  to  wait  until  the  density  drops  near  to 
or  below  the  optimal  density.  They  are  then  free  to  disperse  again  when  the  highway  is 
operating  in  a  non-congested  state. 

Figures  1 1  and  12  represent  similar  cases  to  Figures  6  and  10,  but  with  the  modified  Drake’s 
model  replaced  by  the  modified  Greenshields’  model.  The  optimal  density  of  the  modified 

k  ■ 

Greenshields’  model  is  given  by  k0  =  — .  Since  Ajam  is  set  to  be  200  veh/km  the  optimal 

density  is  therefore  equal  to  100  veh/km.  Similar  results  are  obtained  with  the  modified 
Greenshields’  model.  When  the  initial  platoon  has  a  density  less  than  the  optimal  density, 
dispersion  occurs  throughout  the  analysis  (Figure  11).  If  the  platoon  is  initially  congested 
with  a  density  value  above  the  optimal,  dispersion  is  limited  (Figure  12).  Only  when  the 
density  of  the  platoon  drops  back  to  near  or  below  the  optimal  density  can  vehicles  overtake 
easily  and  the  dispersion  behavior  becomes  clear. 
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4.4  Infinite  number  of  driver  classes 


We  might  wonder  what  physical  meaning  the  small  staircases  carry  in  the  9-class  model.  By 
computing  with  different  number  of  classes  we  have  observed  that  the  number  of  staircases  is 
always  equal  to  the  number  of  classes  (this  corresponds  to  the  different  wave  speeds  of  the 
different  combinations  of  classes),  and  the  strength  of  those  staircases  decreases  with  the 
number  of  classes. 

We  could  thus  consider  an  asymptotic  case  when  the  number  of  classes  goes  to  infinity.  The 
model  (3)  then  becomes  again  a  scalar  equation  but  with  one  more  independent  variable  v, 
corresponding  to  the  distribution  of  driver  classes.  It  reads 

dk(x,t,v)  |  dq(x,t,v)_Q 
dt  dx 


with  the  numerical  flux  given  by 

q(x,t,v)  =  vk(x,t,v)exp[-(k(x,t,v)/ k0)2 /2^j,  for  modified  Drake’s  form  (38) 
and 

q(x,t,v)  =  vk(x,t,v)(l-k(x,t,v)/ kjam),  for  modified  Greenshields’  form.  (39) 

Other  forms  of  traffic  stream  models  can  be  considered  in  similar  fashion.  The  boundary 
conditions  are  now  set  as  a  function  of  the  class  variables  v  e  [vmin,vmax].  We  remark  that 
this  continuous  model  has  some  similarity  with  the  kinetic  models,  however  no  relaxation  is 
involved  and  this  can  be  considered  as  a  relaxed,  equilibrium  model.  The  main  difference  is 
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that,  while  the  conventional  kinetic  models  consider  a  distribution  of  non-equilibrium  speed 
around  an  equilibrium  value,  our  model  assumes  a  continuous  distribution  of  equilibrium 
speed. 

As  example  calculations,  we  assume  a  modified  Drake’s  form  of  traffic  stream  model,  as  in 
equation  (38),  with  ko  =  50  veh/km,  and  an  initial  platoon  of  maximum  density  40  veh/km  as 
shown  in  Figure  1 .  The  distribution  of  k  everywhere  in  the  platoon,  as  a  function  of  v,  follows 
a  continuous  curve  in  the  shape  of  Figure  13,  left,  which,  when  discretized  using  9  points  in 
v,  and  suitably  scaled,  gives  the  original  Figure  4.  Thus,  the  M-class  model  can  be  considered 
as  a  discretization  of  the  continuous  model  (37)  in  the  v  variable.  To  demonstrate  that  the 
solutions  from  the  M-class  model  converge  to  those  of  the  continuous  model  (37),  we  plot  the 
9-class,  21-class  and  41-class  density  versus  distance  graphs  at  t  =  0.015  hours,  in  Figure  13, 
right,  using  WENO  with  1600  points,  which  gives  numerically  convergent  solutions.  We  can 
clearly  see  that  the  solutions  converge  to  smooth  curves  without  staircases  when  the  number 
Mof  classes  increases.  There  is  not  much  noticeable  difference  when  M  increases  beyond  41. 
In  Figure  14,  we  plot  the  density  versus  distance  for  various  times  on  the  left,  and  the  flow 
versus  time  for  various  spatial  locations  on  the  right,  for  the  continuous  model  (37) 
demonstrated  by  the  M  =  41  class  model.  It  is  interesting  to  note  that  a  nice  platoon 
dispersion  behavior  is  observed  for  this  continuous  model. 

5.  DISCUSSIONS  AND  CONCLUSIONS 


The  study  of  traffic  flow  using  a  macroscopic  approach  often  involves  a  single  conservation 
law,  or  a  system  of  conservation  laws  that  are,  in  general,  of  hyperbolic  type.  The  scalar 
LWR  model  and  those  higher-order  continuum  models  proposed  so  far  contain  hyperbolic 
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partial  differential  equations.  Care  must  be  taken  in  solving  these  PDEs  using  numerical 
methods  due  to  the  existence  of  singularities  and  multiple  solutions.  Fortunately  much  work 
on  numerical  methods  for  hyperbolic  PDEs  has  been  carried  out  in  the  field  of  Computational 
Fluid  Mechanics  (CFD),  which  can  also  be  applied  to  traffic  flow  problems.  In  this  paper  we 
applied  one  of  the  state-of-the-art  methods  called  weighted  essentially  non-oscillatory 
(WENO)  scheme  to  obtain  solutions  for  a  recently  proposed  MCLWR  model  by  Wong  and 
Wong  (2001).  The  results  of  a  series  of  numerical  tests  are  encouraging  and  interesting. 
Several  conclusions  can  be  drawn.  First,  for  the  2-class  model,  the  derived  linearized  wave 
speeds  give  reasonable  predictions  of  the  actual  nonlinear  wave  speeds.  Analysis  using 
linearization  on  the  MCLWR  model  also  demonstrated  that  for  this  2-class  case  the  class- 
characterized  waves  never  travel  faster  than  the  fastest  vehicle  in  the  traffic  stream. 

Second,  the  fifth-order  WENO  scheme  has  been  implemented  to  solve  the  MCLWR  model 
and  it  is  more  efficient  than  the  first-order  Lax-Friedrichs  scheme  and  the  first-order 
Godunov  scheme.  The  high-order  WENO  needs  fewer  grid  points  than  the  first-order 
methods  to  obtain  solutions  of  the  same  accuracy.  The  reduction  factor  is  around  64  for  the 
Lax-Friedrichs  scheme  and  8  for  the  Godunov  scheme. 

Third,  from  the  convergence  study  of  the  numerical  methods  used  in  this  paper,  for  first  order 
schemes,  it  might  be  premature  to  accept  that  a  solution  is  converged,  unless  a  very  refined 
mesh  (e.g.  25600  points  for  the  Lax-Friedrichs  scheme)  is  used.  Convergent  solutions  of  the 
evolution  of  initial  platoons  of  vehicles  show  that  the  MCLWR  model  produces  dispersed 
platoons  with  staircase-like  steps.  We  found  that  the  number  of  steps  is  equal  to  the  number 
of  classes  in  the  traffic  stream.  Linearized  analysis  of  the  2-class  model  in  Section  2  shows 
that  there  exist  two  class-characterized  waves,  which  can  also  be  generalized  to  the  M-class 
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model  that  there  would  be  M  different  class-characterized  waves  traveling  in  the  traffic 
stream.  The  speed  of  each  wave  is  characterized  by  the  class-specific  parameters  in  the 
MCLWR  model.  It  is  believed  that  the  formation  of  staircase  is  caused  by  these  class- 
characterized  waves  in  the  traffic  stream.  For  the  2-class  test  case,  each  staircase  is  composed 
of  one  class  of  driver  only.  However,  this  exclusivity  property  does  not  generally  apply  for 
the  cases  with  greater  number  of  classes. 

Finally,  we  have  extended  the  MCLWR  model  to  include  a  continuous  equilibrium  speed 
distribution.  Thus  the  M-class  model  can  be  considered  as  a  discretization  of  this  continuous 
model.  It  has  been  shown,  by  increasing  M,  that  the  continuous  model  can  predict  platoon 
dispersion  behavior  without  the  staircases  observed  in  the  discrete  M-class  model.  This  is 
quite  realistic  as  the  actual  equilibrium  speed  distribution  might  be  expected  to  be  continuous 
in  general.  The  actual  distribution  function  has  yet  to  be  determined  from  field  data;  however 
the  underlying  philosophy  of  the  MCLWR  model  will  not  change. 
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Figure  1 


Initial  density  platoon 
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Figure  3 


The  wave  speeds  for  the  2-class  case. 
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Figure  4 


Distribution  of  drivers  in  the  platoon. 
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Figure  5  Density  versus  distance  at  t  =  0.015  hour.  Top  left:  Comparison  between  first 
order  Lax-Friedrichs  with  6400  points  (solid  line)  and  WENO  with  100  points 
(circles);  Top  right:  Convergence  of  WENO  with  100  points  (dash-dot  line), 
400  points  (dashed  line)  and  1600  points  (solid  line);  Bottom  left: 
Convergence  of  first  order  Lax-Friedrichs  with  400  points  (dash-double  dots 
line),  1600  points  (dash-dot  line),  6400  points  (dashed  line)  and  25600  points 
(solid  line);  Bottom  right:  Comparison  between  first  order  Lax-Friedrichs  with 
25600  points  (solid  line)  and  WENO  with  400  points  (circles). 


Figure  6 


Case  1.  Left:  total  density  as  a  function  of  spatial  location,  for  t  =  0,  0.005,  0.010,  ...  ,  0.025  hours;  Right:  the  flow  as  a  function 
of  time  t  at  x  =  0.2,  0.4,  0.6,  ...  ,2.0  km. 
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Figure  8 


First  order  Godunov  solver.  Top:  convergence  with  mesh  refinements;  bottom: 
density  and  flow  evolution. 
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Figure  1 1 


Case  3.  Left:  total  density  as  a  function  of  spatial  location,  for  t  = 
of  time  t  at  x  =  0.2,  0.4,  0.6,  ...  ,2.0  km. 
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Figure  12  Case  4.  Left:  total  density  as  a  function  of  spatial  location,  for  t  =  0,  0.005,  0.010,  ...  ,  0.025  hours;  Right:  the  flow  as  a  function 
of  time  t  atx  =  0.2,  0.4,  0.6,  ...  ,2.0  km. 
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Figure  14  Continuous  in  v  model  for  M  =  41  classes.  Left:  total  density  as  a  function  of  spatial  location,  for  t  =  0,  0.005, 
0.010,  ...  ,  0.025  hours;  Right:  the  flow  as  a  function  of  time  t  atx  =  0.2,  0.4,  0.6,  ...  ,2.0  km. 
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Table  1  CUP  times  for  the  first  order  Godunov  scheme  and  the  5-th  WENO  scheme 
at  the  same  level  of  accuracy,  (on  a  SunBlade  1000  workstation) 


First  order  Godunov  (nx  =  1600) 

5-th  order  WENO  (nx  =  200) 

t  =  0.1 

20  sec 

4  sec 

t  =  0.2 

39  sec 

8  sec 

First  order  Godunov  (nx  =  6400) 

5-th  order  WENO  (nx  =  800) 

t  =  0.1 

376  sec 

71  sec 

t  =  0.2 

722  sec 

138  sec 

